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We describe in detail full numerical and perturbative techniques to compute the gravitational 
radiation from intermediate-mass-ratio black-hole-binary inspirals and mergers. We perform a series 
of full numerical simulations of nonspinning black holes with mass ratios q — 1/10 and q = 1/15 
from different initial separations and for different finite-difference resolutions. In order to perform 
those full numerical runs, we adapt the gauge of the moving punctures approach with a variable 
damping term for the shift. We also derive an extrapolation (to infinite radius) formula for the 
waveform extracted at finite radius. For the perturbative evolutions we use the full numerical 
tracks, transformed into the Schwarzschild gauge, in the source terms of the Regge-Wheller-Zerilli 
Schwarzschild perturbations formalism. We then extend this perturbative formalism to take into 
account small intrinsic spins of the large black hole, and validate it by computing the quasinormal 
mode frequencies, where we find good agreement for spins |o/M| < 0.3. Including the final spins 
improves the overlap functions when comparing full numerical and perturbative waveforms, reaching 
99.5% for the leading (I, m) = (2,2) and (3,3) modes, and 98.3% for the nonleading (2,1) mode in 
the q — 1/10 case, which includes 8 orbits before merger. For the q = 1/15 case, we obtain overlaps 
near 99.7% for all three modes. We discuss the modeling of the full inspiral and merger based on a 
combined matching of post-Newtonian, full numerical, and geodesic trajectories. 
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I. INTRODUCTION 

There is strong indirect evidence for the existence of 
black holes (BHs) of a few solar masses (M@) residing 
in galaxies and for supermassive BHs (SMBHs), with 
masses 1O 5 M0-1O 1O M0 in the central cores of active 
galaxies. These BHs can form binaries and the merg- 
ers of black-hole binaries (BHBs) are expected to be the 
strongest sources of gravitational radiation and the most 
energetic event in the Universe. The current generation 
of ground-based interferometric gravitational wave detec- 
tors, such as LIGO, VIRGO, and GEO, are most sensi- 
tive to BHB mergers with total masses of a few tens to 
hundreds of solar masses, while the space-based LISA 
detector will be sensitive to mergers of BHBs with a few 
million solar masses. 

The existence of intermediate-mass BHs (IMBH) , from 
a few hundred to tens of thousand of solar masses, is 
still uncertain. If they exist, then these IMBH can 
form binaries with solar-mass-sized objects, leading to 
compact-object mergers with mass ratios in the range 
0.001 < q = mxlm-x < 0.1, which could be detected 
by advanced LIGO. The detection of gravitational waves 
from these encounters, as well as the correct modeling of 
the waveform as a function of the BHBs physical param- 
eters, would allow us to estimate the population of such 
objects in the Universe. Likewise, encounters of IMBH 
with SMBHs in the centers of a galaxies would lead to 
mergers with mass ratios in the range 0.001 < q < 0.1, 
detectable by LISA. On the other hand, theoretical N- 
body simulations [T], assuming direct cosmological colli- 



sions of galaxies with central SMBHs, set the most likely 
SMBH binary mass ratios in the range 0.01 < q < 0.1. 

In Refs. [H [3] the prospects of detecting IMBH bi- 
nary (IMBHB) inspirals with advanced LIGO was dis- 
cussed, and in Ref. [4] it was shown that intermediate- 
mass-ratio (IMR) inspirals of IMBHs plunging into su- 
permassive BHs would be relevant to LISA, while IMR 
mergers of IMBHs with stellar objects can be detected by 
LIGO/VIRGO. In both cases the accuracy of the post- 
Newtonian (PN) approach (which was used to model the 
gravitational radiation) was questioned and the need for 
more accurate waveforms was stressed. 

After the 2005 breakthroughs in numerical relativity 
[5H7], simulations of BHBs became routine. The explo- 
ration of generic binaries [5] led to the discovery of large 
recoils acquired by the remnant BH. While long term 
generic BHB evolutions are possible, including the last 
few tens of orbits [HI HH] , two very interesting corners of 
the intrinsic parameter space of the BHBs remain largely 
unexplored: maximally spinning binaries and the small 
mass ratio limit. 

In a previous letter [TT] we introduced a new technique 
that makes use of nonlinear numerical trajectories and 
efficient perturbative evolutions to compute waveforms 
at large radii for the leading and nonleading modes. As a 
proof-of-concept, we computed waveforms for a relatively 
close binary with q — 1/10. In this paper we will describe 
these techniques in detail, extend them to slowly spinning 
black holes, and reach smaller mass ratios, to the q = 
1/15 case, with full numerical simulations. 

The paper is organized as follows. In Sec. [H] we de- 
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scribe the full numerical techniques employed in the evo- 
lution of BHBs. Those are based in the moving punc- 
ture approach [5] [B] with a gauge choice that allows a 
spatial and time variation of the gamma-driver parame- 
ter rj(x a ,t). We describe the results of the simulations 
for two different mass ratios q = 1/10,1/15 and two 
different initial separations leading to evolutions with 
BHs performing between 4 and 8 orbits prior to merger, 
the latter representing the longest waveform published 
so far in the small q regime. The gauge has also been 
shown to work for evolutions of a nonspinning q = 1/100 
BHB [12]. In Sec. Ill we describe the perturbative tech- 
niques used to evolve a particle around a massive black 
hole. We extend the Regge-Wheeler-Zerilli (RWZ) tech- 
niques to include, perturbatively, a term linear in the 
spin of the larger black hole. This takes the form of 
second-order perturbations and adds a source term to 
the usual Schwarzschild perturbations (SRWZ). We also 
study the asymptotic behavior of the perturbative solu- 
tions for large r and come up with a practical way of 
correcting finite observer location effects perturbatively 
on the numerical waveforms. In Sec. IIVI we describe the 
results of comparing full numerical waveforms with per- 
turbative ones that use the full numerical tracks for the 
particle motion. We compute matching overlaps for the 
leading modes {I, m) = (2,2); (2,1); (3,3). We verify 
the scaling of the waveform amplitudes with the reduced 
mass n for the mass ratios q= 1/10, 1/15. We also quan- 
tify the effects of adding the spin of the final black hole 
into the perturbative integrations. In Sec. |V]we discuss 
the properties of the full numerical trajectories in the two 
cases studied q = 1/10,1/15 that can be generalized to 
smaller mass ratios and hence can help in providing a 
modeling for the tracks used in the perturbative integra- 
tion, in particular, the final "universal plunge" and the 
use of resummed PN trajectories for the stages prior to 
the full numerical simulation. Finally in the Appendix [A] 
we give further evidence of the accuracy and validity of 
the SRWZ formalism here developed by computing the 
quasinormal modes (QNM) and comparing them with 
the exact Kerr black-hole modes for different values of 
the spin parameter. 



II. NUMERICAL RELATIVITY TECHNIQUES 

To compute the numerical initial data, we use the 
puncture approach [13] along with the TwoPUNC- 
tures [T3] thorn. In this approach the 3-metric on the 
initial slice has the form y a f, = (ipBL + u) 4 8ab, where 
ipBL is the Brill-Lindquist conformal factor, S a b is the 
Euclidean metric, and u is (at least) C 2 on the punc- 
tures. The Brill-Lindquist conformal factor is given by 
tpBL = 1+X)j=i m i/(2\r — fj|), where n is the total num- 
ber of 'punctures', mf is the mass parameter of puncture 
i (rnf is not the horizon mass associated with puncture i), 
and fi is the coordinate location of puncture i. We evolve 
these black-hole-binary data-sets using the LazEv [TS] 



implementation of the moving puncture approach [S] |B] 
with the conformal factor W = y/x — exp(— 20) sug- 
gested by [TB] For the runs presented here we use cen- 
tered, eighth-order finite differencing in space [17] and 
an RK4 time integrator. (Note that we do not upwind 
the advection terms.) 

We use the Carpet [18] mesh refinement driver to pro- 
vide a "moving boxes" style of mesh refinement. In this 
approach refined grids of fixed size are arranged about 
the coordinate centers of both holes. The Carpet code 
then moves these fine grids about the computational do- 
main by following the trajectories of the two black holes. 

We use AHFinderDirect [TH] to locate apparent 
horizons. We measure the magnitude of the horizon 
spin using the Isolated Horizon algorithm detailed in 20 . 
This algorithm is based on finding an approximate rota- 
tional Killing vector (i.e. an approximate rotational sym- 
metry) on the horizon ip a . Given this approximate Killing 
vector tp a , the spin magnitude is 



S l<p] = 



~ / (<p a R b K ab )d 2 V, 
<«r J AH 



(1) 



where K ab is the extrinsic curvature of the 3D-slice, d 2 V 
is the natural volume element intrinsic to the horizon, 
and R a is the outward pointing unit vector normal to 
the horizon on the 3D-slice. We measure the direction of 
the spin by finding the coordinate line joining the poles 
of this Killing vector field using the technique introduced 
in [21] . Our algorithm for finding the poles of the Killing 
vector field has an accuracy of ~ 2° (see [2T] for details) . 
Note that once we have the horizon spin, we can calculate 
the horizon mass via the Christodoulou formula 



,H 



m? r + ^/(4mf rr ), 



(2) 



where m- ur = \J Aj (167r) and A is the surface area of the 
horizon. We measure radiated energy, linear momentum, 
and angular momentum, in terms of using the for- 
mulae provided in Refs. [22], [23] . However, rather than 
using the full ip4, we decompose it into £ and m modes 
and solve for the radiated linear momentum, dropping 
terms with I > 5. The formulae in Refs. [22] [23] are 
valid at r = oo. Typically, we would extract the ra- 
diated energy-momentum at finite radius and extrapo- 
late to r = oo. However, for the smaller mass ratios ex- 
amined here, noise in the waveform introduces spurious 
effects that make these extrapolations inaccurate. We 
therefore use the average of these quantities extracted at 
radii r = 70, 80, 90, 100 and use the difference between 
these quantities at different radii as a measure of the er- 
ror. We found that extrapolating the waveform itself to 
r = oo introduced phase errors due to uncertainties in the 
areal radius of the observers, as well as numerical noise. 
Thus when comparing perturbative to numerical wave- 
forms, we use the waveform extracted at r — 100A/. In 
Sec. |IIIB 7| we provide an alternative method of extrapo- 
lation of waveforms based on a perturbative propagation 
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of the asymptotic form of 7^4 at large distances from the 
sources leading to the following simple expression 



lim [r^ m (r,t)] 



-o(Rot) 



(£-!)(£ + 2) 



dt^ e 4 m (r,t) 



(3) 



where robs is the approximate areal radius of the sphere 
-Robs = const [Add a factor (1/2 — M/r) multiplying the 
square bracket to correct for a difference in normalization 
between the Psikadelia and Kinnersley tetrads at large 
distances.] We have found that this formula gives reliable 
extrapolations for i?obs ^ 100M. 



where we chose Rq = 1.31. The above gauge condition 
is inspired by, but differs from Ref. [25] between the BHs 
and in the outer zones when a / 1 and b =/= 2. Once 
the conformal factor settles down to its asymptotic tp — 
C/y/r + 0(1) form near the puncture, rj will have the 
form rj = (R /C 2 )(l + b(r/C 2 ) a ) near the puncture and 
i] = R r b ~ 2 M I (aM) b as r — > 00. In practice we used 
a = 2 and 6 = 2, which reduces 77 by a factor of 4 at 
infinity when compared to the original version of this 
gauge proposed by [25] . We note that if we set b = 1 then 
7] will have a 1/r falloff at r = 00 as suggested by [25] . 
Our tests indicate that the choices (a = 2, b = 1) and 
(a = 1,6 = 1) lead to more noise in the waveform than 
(a = 2,6 = 2). 



A. Gauge 

We obtain accurate, convergent waveforms and horizon 
parameters by evolving this system in conjunction with 
a modified 1+log lapse and a modified Gamma-driver 
shift condition [3 [24], and an initial lapse a(t = 0) = 
2/(1 + ip% L ). The lapse and shift are evolved with 

(d t -l3 l di)a = -2aK, (4a) 
8 t /3 a = (3/4)f a -r,(x a ,t)(3 a , (4b) 

where different functional dependences for n(x a ,t) have 
been proposed in [15, 25 29 . Here we use a modification 
of the form proposed in [25], 



n(x a ,t) = R Q 



v / d l Wd J Wfi 
(1- W a ) b 



(5) 



B. Simulations and results 



In order to obtain low-eccentricity initial data parame- 
ters, we started with quasicircular post-Newtonian initial 
data parameters for the momenta and particle positions. 
We then evolved for 1-2 orbits, and used the procedure 
detailed in 30 to obtain lower eccentricity parameters. 
In practice we performed between 3 and 4 iterations of 
the above procedure. In Table [I] we show the initial data 
parameters, horizon masses and mass ratio, and initial 
orbital eccentricities for the three configurations consid- 
ered here. 



TABLE I: Initial data parameters. The punctures are located on the x-axis at positions x\ and X2, with puncture mass 
parameters (not horizon masses) mi and m2, and momentum ±p. In all cases, the punctures have zero spin. Configuration 
ql0r7.3PN is based on the original PN parameters, prior to any eccentricity removal iteration. The lower part of the table 
shows the horizon masses mn, and mH 2 , the mass ratio q, the ADM mass, and the eccentricity e. 



Config 


Xl 


X2 


Px 


Py 


m 1 


ni2 


gl0r8.4 


7.633129 


-0.7531758 


-0.000168519 


0.0366988 


0.08523727 


0.90739686 


gl0r7.3 


6.604383 


-0.6715184 


-0.000219713 


0.0410386 


0.08438951 


0.90703855 


ql0r7.3PN 


6.604383 


-0.6715184 


-0.000326708 


0.0404057 


0.08438951 


0.90703855 


gl5r7.3 


6.806173 


-0.4438775 


-0.000160518 


0.0290721 


0.05756623 


0.93622418 


Config 




m H2 


q 


A^ADM 


e 




gl0r8.4 


0.091289 


0.912545 


0.10004 


1.0000428 


0.0004 




?10r7.3 


0.091378 


0.913010 


0.10008 


1.00025882 


0.0017 




gl0r7.3fW 


0.091329 


0.912990 


0.10003 


1.00000000 


0.008 




g!5r7.3 


0.062536 


0.940421 


0.06650 


1.00005083 


< 0.0015 





In all the simulations presented here, the outer bound- three resolutions, with a global refinement factor of 1.2 
aries were placed at 400Af. We performed runs with between resolutions. For the q — 1/10 runs, we used 
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TABLE II: Remnant horizon parameters and radiated energy-momentum 



Config E rad Jrad M H - Madm Sadm - S H a Kick km s 

gl0r8.4 0.00446 ± 0.0001 0.0517 ± 0.001 0.00046 ± 0.00003 0.05028 ± 0.00001 0.25986 ± 0.00001 59.4 ±3.0 

gl0r7.3 0.00400 ± 0.00001 0.0386 ± 0.003 0.00415 ± 0.00001 0.04028 ± 0.00001 0.26034 ± 0.00001 65.8 ± 2.0 

gl5r7.3 0.00216 ± 0.00001 0.0235 ± 0.0004 0.00225 ± 0.00001 0.02289 ± 0.0004 0.18872 ± 0.00001 33.5 ±2.1 



11 levels of refinement around the smaller BH, with a 
central resolution of h = M/307.2 for the coarsest runs, 
while for the q — 1/15 run we used 12 levels of refinement, 
with a central resolution of M/614.4. In Table[ll]we show 
the radiated energy-momentum and remnant BH param- 
eters for these configurations. In the figures and tables 
below we refer to the different resolution runs using the 
gridspacing on the coarsest grid relative to ho = 10/3M. 

FIG. 1: The puncture separation as a function of time for 
three q = 1/10 simulations. The solid curve shows a high- 
eccentricity simulation obtained from PN quasicircular pa- 
rameters (with particle limit corrections); the dotted curve 
shows results from a simulations with similar initial separation 
after a few iterations to reduce eccentricity; the dot-dashed 
curve shows an even further separated binary with still smaller 
eccentricity. Note that the initial jump in the orbit does not 
appear to be a strong function of the eccentricity or initial 
separation. 
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In Fig. [T] we show the orbital separation as a function of 
time for the <?10r8.4 and ql0r7.3 configurations, as well as 
a high-eccentricity configuration obtained by directly us- 
ing PN parameters in the initial data (ql0r7.3PN) that 
we used for the proof-of-concept in Ref. [TT]. Note that 
the highly eccentric ql0r7.3PN binary merges sooner 
than the lower eccentricity ql0r7.3. From the plot we 
can also see that the initial jump in the orbit is not a 
function of either initial separation or eccentricity. In 
Fig. [2] we compare the orbital separation for the t7l0>7.3 
and c7l5r7.3 configurations. From the plot it is clear that 
the initial jump in the orbit is not a strong function of 
mass ratio either. This indicates that the initial jump 



FIG. 2: The magnitude of the puncture separation (|a?i — a?2 1 ) 
as a function of time for a q = 1/10 and q = 1/15 binary 
at similar initial separations. Note that the initial jump in 
the orbit appears to be independent of q. Also note that the 
q = 1/15 run inspirals more slowly. 
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will become more problematic as the mass ratio is re- 
duced (and hence the inspiral becomes weaker). We also 
observe that, quite independent of the initial separation 
and the initial eccentricity, the track displays a univer- 
sal behavior during the final plunge. This confirms that 
the tracks are gravitational radiation driven and we are 
numerically resolving this radiation accurately. 

In Fig.[3]we show the orbital trajectories of the gl0r7.3 
and gl5r7.3 configuration. In the plot the curves have 
been rotated to maximize the overlap during the plunge. 
From the plot we see a "universal" plunge behavior at 
small separations with distinctly different orbital dynam- 
ics at larger separations. As expected, the small mass 
ratio binary merges more slowly. In Fig. [4] we show the 
real part of the (£ = 2, to = 2) mode of ip4 for the t7l0>7.3 
and (?15r7.3 configurations. Here the we rescaled 1P4 for 
ql5r7.3 by a factor of 1.5. Note that the good overlap of 
the rescaled ip4 indicates that the amplitude of ip4 scales 
with q (before the different orbital dynamics of q = 1/10 
and q — 1/15 cause the ql0r7.3 to merge sooner). In 
Fig. [5] we show the convergence of the ylO>7.3PA?~ con- 
figuration for three resolutions. Note that in this plot, 
the low resolution actually corresponds to a grid-spacing 
1.2 times larger than the low resolutions for the other 
configurations. From the plot we can see that at later 



5 



FIG. 3: An (xy) projection of the puncture separation (xi — 
X2) for a q = 1/10 and q — 1/15 binary at similar initial 
separations. The trajectories have been rotated so that they 
overlap during the plunge and merger. Note the "universal" 
plunge trajectory. 
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FIG. 4: The real part of the (£ = 2, m = 2) mode of i/>4 
for a q — 1/10 and q = 1/15 binaries starting at similar 
separations. The waveform from the q — 1/15 binary was 
rescaled by a factor of 1.5 (15/10). 
0.0001 

V 4 (q10r7.3) 
(15/10) *v>„(q15r7.3) 



5e-05 - 



S 



-5e-05 - 



-0.0001 





200 



400 
t/M 



600 



800 



FIG. 5: The convergence of the phase and amplitude of 
the (£ = 2,m = 2) mode of ip A for the qWr7.3PN config- 
uration. Note that here the three resolutions consist of a 
low resolution with grid-spacing 1.2 times larger than the 
low resolution runs for gl0r7.3, ql0r8.4, gl5r7.3 configura- 
tions. Eighth-order convergence implies ^4(1.2/10) — ip4,(ho) = 
4.29982(i/>4(/io) — il>4(ho/1.2)), while fourth-order convergence 
implies 4> 4 (1.2h ) - ipi{h ) = 2.0736(^ 4 (^o) - ipi{h /l.2)). 
Initially, the error in ip4 is very small and dominated by grid 
noise. Eighth-order convergence in the amplitude is apparent 
beginning at t = 320M, while eighth-order convergence in the 
phase becomes apparent at t = 420M. The dashed vertical 
line shows the time when the wave frequency is Mu = 0.2. 
The phase error at this frequency is 5(j> < 0.2 rad. 
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time the convergence is eigth-order. The earlier time 
fourth-order convergence is due to finite-difference and 
interpolation errors in the extraction routines. At later 
times, the phase error dominates the errors in the wave- 
form, and this error converges to eighth-order. Finally, 
in Fig. [6] we show the phase of the waveform for qlbr7.3 
for three resolutions. The phase errors near the plunge 



are reported in Table III 



III. PERTURB ATIVE TECHNIQUES 

In this section we describe in some detail the use of 
perturbative techniques to produce BHB waveforms from 
a small mass ratio system. We summarize the key for- 
mulae used (for more details see, for instance, [31]). and 
extend the formalism to add the spin of the large black 
hole as a second-order perturbation, coupling it to the 
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FIG. 6: (Top) The phase of (t = 2, m = 2) mode of 4> 4 
for a q = 1/15 BHB for three resolutions. Note that the 
phase error only converges to fourth-order and that the high- 
est resolution is refined by a factor of 1.2 2 rather than 1.2 
with respect to the medium resolution. (Bottom) A con- 
vergence plot showing the initial (better than) fourth-order 
convergence of the waveform. Note here that the differences 
V>4(l-2ft ) - ipi{h Q ) = 1.39895(^ 4 (/io) - ^ 4 {h /1.2 2 )) if the 
waveform is fourth-order convergent. 
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radiative first-order perturbations. We neglect quadratic 
terms in the radiative modes of the order 0{q 2 ). The 
resulting equations are still of the Regge- Wheeler and 
Zerilli form (we are still doing perturbations around a 
Schwarzschild background), but they now include ex- 
tended source terms with linear dependence on the spin 
in addition to the local (Dirac's deltas) source terms 
already present in the first-order formalism. We plug 
into these latter terms the full numerical trajectories 
(hence indirectly also adding a spin dependence). We 
denote the resulting formalism as Spin-Regge-Wheeler- 
Zcrilli (SRWZ). 

A. Metric perturbations and particle's orbit 

1. Spin as a perturbation 

We consider the Kerr metric up to 0(a 1 ). Here a de- 
notes the spin of the black hole which has the dimension 
of mass. In the usual Boyer-Lindquist coordinates, this 
is given by 

2 r - 2 M 2 Ma sin 2 6d(j) dt r 2 



+r 2 d6 2 + r 2 sin 2 9d(j) 2 + 0(a 2 ) . (6) 

In the above metric, the terms which depend on a 
are treated as the perturbation in the background 
Schwarzschild spacetime. 

9»» = glf + /^ spin) ■ (7) 

For the above metric perturbations, we consider the ten- 
sor harmonics expansion defined using the tensor har- 
monics of [35] . We find that the first-order perturbation, 
0(a 1 ), is related to the I — 1, m = odd parity mode, 
and the coefficient of the tensor harmonics is given by 

^o pin) (^) = Vtt- (8) 

where S — Ma. The other components are zero. 



TABLE III: The phase error in the [l = 2, m = 2) mode of ip 4 
[extracted at R — 100M and extrapolated to oo using Eq. |3|] 
when the waveform frequency is Mui = 0.2 for the medium 
and high-resolution runs. The table shows the predicted phase 
errors extrapolating to infinite resolution and assuming eigth- 
and fourth-order convergence. 



Config 




Eigth-order 


Fourth-order 


gl0r8.4 (h = 
gl0r8.4 (h = 
gl0r8.4 (h = 
gl5r7.25 (h = 


ho) 

h /1.2, pred) 
ho/1.2 2 , pred) 
= ho/1.2 2 ) 


0.205133 
0.0477073 
0.0110952 
0.1406 


0.630496 
0.304058 
0.146633 
0.762 



2. Second-order formulation 

In the following, we treat spin-radiation couplings in 
the second-order perturbation. Therefore, we consider 
the Einstein equation in the second perturbative order. 

= 8 it (tW+T< 2 ,)) =8nT^, (9) 

According to [33] , and the fact that we use the Numerical 
Relativity (NR) trajectory, we do not separate the first 
and second-order energy-momentum tensor of the parti- 
cle. And the second-order metric perturbation, /j( 2 > wavc ) 
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is created by the spin, /i (1 ^ spin) -radiation, h^' wavc ^ cou- 
plings. In this case, we may solve 



gM^wavc)] = g7rT ^ ) (1Q) 
G (l)^(2,wavc)] = _ G (2)^(l,wavc)^(l, S pin)^ ^ 



up to 0(a 1 ), where we ignore the square of the first-order 
wave functions. 



As discussed below, we solve Eqs. (10) and (11) for 



the even parity perturbation of the Regge-Wheeler-Zcrilli 
formalism in the following form. 

G$[/i (1 ' wavo) + /i( 2 > wavc )] + G$[h {1 '™ ve \ h {1 ^] 
= GW[h (wavc) ] + G$[/i (wavo) , h^ spin) ] 
= SnT^, (12) 

where M wavc ) = h^' wave ^> + /j(2,wave)_ On the other hand, 



for the odd parity perturbation, Eqs. (10) and (11) are 
solved for each perturbative order. 

Here we consider intermediate mass ratio binaries. As 
discussed in we can introduce some second-order 

effects that arise purely from the particle's first-order 
perturbation, if we treat the particle as a reduced mass 
fi = m\mil {m\ + m-z) orbiting around a black hole with 
the total mass M — mi + mi . 



3. Orbit for inspired 

First, we should note that the coordinates used in NR 
simulations are chosen to produce stable evolutions and 
correspond, initially, to isotropic coordinates. Perturba- 
tive calculations, on the other hand, regularly make use of 
the standard Schwarzschild coordinates. The easiest way 
to relate the two is to translate the numerical tracks into 
the Schwarzschild coordinates. This can be achieved by 
considering the late-time numerical coordinates that cor- 
respond to radial isotropic "trumpet" stationary 1 + log 
slices of the Schwarzschild spacetime [35] . We obtain the 
explicit time and radial coordinate transformations fol- 
lowing the procedure detailed in Ref. (36) . 

Thus, we consider the NR trajectory as an orbit pro- 
jected on the Schwarzschild background. Therefore, we 
calculate the particle's energy, angular momentum etc. 
by using the Schwarzschild metric. Here, since we have 
only the three velocity v l (t) from the data of the NR 
trajectory, the time component of the four velocity is 
derived by assuming the "instantaneous" Schwarzschild 
geodesic approximation. 

In this approximation, the energy and angular momen- 
tum are given by. 



E 



2M\ , 

1 u . 

R J 



R 2 v* 



(13) 
(14) 



where u M = dx^/dr is the four velocity, R = R(t) denotes 
the orbital radius, and we are considering the equatorial 



orbit (9o = tt/2). To evaluate U(t) = it*, we use 
g^u^u" = -1 



{U{t)f 



2M 

W)J 



Wi 



2M 
W) 



-(R(t)f($(t)f 



(15) 



Here, R = u r ju l = dR/dt and $ = ju l = d<S>/dt are 
the three velocity of the particle. 

We note that the energy E derived from the above 
U(t) does not decrease monotonically, and also in the 
end of the orbital evolution, we can not calculate U(t) 
appropriately by using Eq. (15), because U(t) — > oo or 
becomes complex. U(t) — > oo is, in practice, not inconsis- 
tent because U(t) - (1 - 2M/R(t)y 1 for Schwarzschild 
geodesies. 

Therefore, we fix the energy at some orbital radius (or 
time t = t m ) as 



2M 

R(t m ) 



U(t m ) , 



(16) 



and use the following expression to obtain U(t) for 
smaller radii. [This may give the innermost stable cir- 
cular orbit (ISCO) radius.] 



U(t) 



E„ 



1 



2M\ 1 



(17) 



At this stage, we still use the three velocity derived from 
the NR trajectory. 

Here we set R(t m )/M = 7.64 for the q = 1/10 case. 
This radius is obtained from the energy minimum evalu- 
ated by Eq. (14). In the q = 1/15 case, we do not have 



such an energy minimum. Therefore, we simply set the 
same radius as for the q = 1/10 case. 



4- Orbit near merger 

There are large differences between the coordinate sys- 
tem used in the NR simulation and the Schwarzschild co- 
ordinates near the horizon. Although the binary merges 
at finite time in the NR simulation, the binary does not 
merge in the Schwarzschild coordinates. Therefore, we 
need to give the orbit near the horizon. 

Here, we assume that the radiation reaction is not im- 
portant near merger after t — tf, and use the geodesic 
orbit on the Schwarzschild spacetime. First, we consider 
the conserved quantities, i.e., the energy and angular mo- 
mentum. 



E, n — E(t m ) — E(tf) . 
Lf = L{t f ) 

= R(t f ) 2 ^>(t f )E n 



2M V 



(18) 
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where E m is the same as the previous section. And then, 
from the above equations, we calculate 



U(t) 



1 



2M 



R(t) 
L f R(t) - 2M 



(19) 



On the other hand, we use a fitting formula for the 
radial trajectory. By using g^u^u 11 = —1, we define an 
effective energy for the radial motion, 



, / 2M Y 



J2 



R(t f Y 



R{tff 



(20) 



and consider E r as a constant after t = tf. The evolution 
of R(t) is derived as 

( 2M \ 



\ 



l 

m 



2M \ 

W)J 



R(ty 



(21) 



From this equation, we can obtain various equations if we 
need, for example, R(t) = (8R(t) / dR{t))R{t) . It is noted 
that we may consider another treatment as discussed in 
Sec.lV] 



In our perturbative code for both q — 1/10 and 1/15 
cases, we set R(tf)/M = 3.0 which is inside the ISCO 
radius. This is because we want to use the NR trajecto- 
ries as long as possible in this paper, and the data of the 
tracks become noisy inside the above orbital radius due 
to the coordinate transformation. 



B. Regge-Wheeler-Zerilli equations with spin 

1. First-order Regge-Wheeler-Zerilli equations 

For the notation of the Regge-Wheeler-Zerilli formal- 
ism [371 EE] , we use [32] and [34] . In the first-order per- 
turbation, i.e., the nonspinning case, we may solve the 
equations, 



dt 2 
-V. 



(even) / (. \ _ r-(even,l) 



(r)*^(t,r) = 5^(t,r) ) (22) 



for the even parity with the Zerilli function ^f^, and 

-^ (odd) (r)*& 1 J(t,r) = S^ d ' 1 >(*,r) 1 (23) 
for the odd parity with the Regge- Wheeler function 

Herer* = r + 2M ln[r/(2M) - 1] is a characteris- 
tic coordinate, and the first-order source terms, sj^ en>1 ^ 
and S^f d are given by 



S 



I'm 



■1) 



(t,r) 



S. 



(odd.l) 
Im 



(*,r) 



167r(r-2M) 2 (r£ 2 +rf-4r + 2M) (1) 10 \ &-k\t - zivi) (1) 

£(£+l)(r£ 2 +r£-2r + 6M)r < m( *' T > ~ ~J?(e + - l)(/+2) £m 



32 7 r(r-2M) 2 v / 2 



16v^?r(r-2M) 

w 

32 vr(r - 2 M) 3 



(r£ 2 + rl - 2 r + 6 M) y/Kjt + 1) ~ lm ' 7 {r£ 2 + rl - 2 r + 6 M)l(l + 1) Or £mK ' ; 
16 tf r(£ 4 r 2 + 2 r 2 £ 3 - 5 r 2 £ 2 + 16 r£ 2 M ~6r 2 £+ 16 r£M + 8 r 2 - 68 rM + 108 M 2 ) (1) 



32 ?r(r- 2 M)r 2 
16\/2 7r(r-2M) 



+ l)£(r£ 2 +r£-2r + 6 M) 2 
d .m ,. s 32V2tt(7 



,(i) -j , azy^ l r-2M) 2 (1) 

{rl 2 + r£-2r + 6 M)£{£ + 1) dr otmK ' y (rf 2 + r£ - 2 r + 6 M )^ + 1) ern 



l6V2inr(r - 2 M) 



,(D (t r) , 1672 7rr(r-2M) ,9 (1) 
W ' j ^(ITT)^ - + 2) dr m 



vWTT)(£- i)(£ + 2) 3* 



a ea(*.r), 



(24) 



where, A^l etc. denote the tensor harmonics coeffi- odd parity wave function ^T^ 1 are related to the Mon- 



(o,l) 



cient of the particle's energy- momentum tensor T^ v . It 
is noted that the even parity wave function ^2 and 



crief 's [5S] and the Cunningham et al. 
a normalization factor, respectively. 



waveforms by 



9 



2. Even parity perturbation with spin where the second-order source term S^ en in the above 

equation is given by 

When we discuss only the second-order Einstein equa- 
tion in Eq. (Ill for the even parity perturbation, the 
Zerilli equation with the 0(a 1 ) spin effect is written as 

- — M 2) it r) + — * (2) (t r) 



Vt en \r)*fj l (t,r) = st en ' 2) (t,r), (25) 



Sir^ (t, r) = 5t v ° n ' 2) (E, S) + 5t ven » 2) (0,5); 

.(cve„,2) r 9 , = rnS / V2t: (-r + 2 M) (-2 r + rl + rl 2 + 12 Af ) (1) 

fm 1 ' j f(f + l)(rP+rf-2r + 6M) V y/t{t+\)r(rP+rl-2r + GM) Mm 



V27rJ-r + 2M) (1) 192 i (— r + 2 M) 7r \/2 (1) 



8»(12M -6r + ^r + 2rf 3 + r£ + 2r£ 2 ) (-r + 2M) ^ ^ | 8il(l + l) fl ^(i) 
r 3 ( r ^2 + r ^ _ 2 r + 6 M) 



-5 r 2 ^ - 28 rM + 12 r£M + 4 f 3 rM + 12 r£ 2 M + 36 M 2 ) (^ - 1) (I + 2) /[(rf 2 + rl - 2 r + 6 Af) 2 r 5 ] 
x (r dth ( ll lm (t, r) + 2 h$_ lm (t, r) - r d r h^_ lm (t, r)) 



v^tt (^ + 2) (r^ 2 + rl - 2 r + 3 M ) (^ - l) 2 (r - 2 M) (1) 

r ^(^_ 1 )^( r ^2 + rf _2r + 6Af) 2 _Q °^ lm ' 



+ ^ 4 / ^+m + lH£-m+l) / _ 2 2 _ 2 

+ ^£+l)(£-l)(£ + 2)\/ (2£ + l)(2£ + 3) \ b{r iM )Vrl+6lrM 6bM 

+3 l 4 r 2 + 4 f 3 rM + r 2 £ 5 + 2 r 2 ^ 3 - 8 r 2 ) - 1) {I + 2) /[(r£ 2 + rl - 2 r + 6 Af) 2 r 5 ] 

x (r ftftiV+im (*> O r + 2 Z^+i™ («, r) - r <U«+i m (*, r)) 



V2tt (I -I )(rl 2 +rl-2r + 3M) (1 + 2 2 (r - 2 M) m , 

-32 V ; n = ^ QA7 +lm (t, r . 26 

r V(£ + 1) {I + 2) {rl 2 + rl-2r + 6 M) 2 ' ' 



SZ en ' 2> {E, S) and S£° D ^{O, S) mean the coupling be- 
tween the black hole's spin and the first-order even and 
odd parity perturbations, respectively. The tensor har- 
monics coefficients of the first-order metric perturbation, 
ffjl^ etc. are written in terms of the first-order Regge- 
Wheelcr and Zerilli functions. 



Here, we introduce the following combined function. 

W*>r) = *&(t,r) + *$(t,r) , (27) 

which is the linear combination of the first- and second- 
order wave functions. This function formally satisfies 
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d 2 T , v a 2 



r*A» ( f > r ) + 1T3*A» (*» r ) " V t Cn) (O**™ (t, r) 



+iSm P { r cn ' 1] (r)^ im (t,r)+iSm P £ (cvcn < 2) (r) ^ 9 tm (t, r) 



= S 



I (l- m ){l + rn [_ (oven,-) (o,!) 
V (2^-1) (2^ + 1)^ W^-lm^.H 

^ t !j (2/" + W 3)' ^ ' +) ^ '> + ^ ' ( 2§ ) 



where iS& ren ' L ' (t, r) denotes the local source term with 
the Dirac's delta function and its derivative. The ex- 
plicit expression and some detailed analysis are given in 
Appendix |A 1| 



3. Odd parity perturbation with spin 



(2) (2) 

Q e ^ and T>\' are calculated by the tensor harmon- 



ics expansion of -G]fJ [/j,( 1 » wave ) , /i( 1 ' spin )]/(87r) from the 

(2) 

second-order Einstein tensor. And using Q\ ' we have 



the relation between the two waveforms \I/ 



*£: z ' 2) as 



(o,2) 
('rn 



and 



In the first perturbative order calculation, we have 
used the Cunningham et al. waveform v E'j>°„ 1 ' ) for the 
odd parity as the Regge- Wheeler function. When we use 
we have some trouble in the source terms of the 



(o,2) 
tin 



perturbed Regge- Wheeler (odd parity) equation. The 
second-order local source term does not vanish at the 
horizon. Therefore, we use the Zerilli waveform / "' 

instead of the Cunningham et al. waveform "f^f 
second perturbative order 



in the 



d 2 



<f' 2) (*,r) + 



if 



Q t 2 -tm \"fJ ' g r «2 

-Vt M \r)^' 2) (*,r) = S^ 2) (*,r) , (29) 

where the second-order source term S^ d ' Z ' 2 ' > is formally 
given as 



c,(odd,Z,2j , 



T 2 y/l (£+1 

I 

8V2TTi(r-2M) 2 



16V2mM(r~2M) p(2) 
r 2 ^i(£ + l)(£-l)(£ + 2) £ml ' 



rv^TI)FW+2j* 



9 ^((/). (30) 



16 \Z2iri r ( 



2*& z ' a, (t 1 r) 



2M) 



(31) 



For the wave equation of 
order source term as 



, we have the second- 
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St**™ it, r) = (E, S) + (0,5); 



-12 



45 



£(£ + !)(£ -1) (£ + 2) ']] (2t-l) (2£ + l) 
V2ir (r-2M) (£-l)(£ + 2) (t + 1) (1) 



(l-m)(l + m) / (r-2M) (*-!)(* + 2) (l+l) fl jA1) 

6 ~a °t JX e.-lm l r > r ) 



r 2 ^{£~l)£(£-2) {£+!) 



45 



^ + l)(£-l)(£ + 2)y (2^+l)(2£ + 3) 
V2tt (r-2M)£(£-l)(£ + 2) (1) 



(£ + m + l)(£-m + l) f , (r-2M)<(M)(<+2) ay(1) 

„4 a t ft £+lra I 1 ' ~J 
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rV(Hl)(f + 2)f(H3) 



,(odd,z,2 )(0 ^ Sm f -48V2nr (r-2M) n(1) 



12i(6r + r* + 7*-14M)(r-2M) (1) 



+ 1) V (^ + i)r 3 

12»(r-2M) (3r-7M) (1) 4 j (-9 r + + rg + 21 M) (r - 2 M ) (1) 

~6 "rfl-Mm ^ r i H „6 ^"Wm ^' H ) > 



(32) 



without any regularization (or modification) of the wave 

function. Here, we note that 5^ d,Z ' 2 ^ (E, 5) for the £ = 
3, m = 2 mode is the time derivative of the second- 
order term in Eq. (76) of [STj . The explicit expression of 
Eq. ( 29 ) is given in Appendix A 2\ We should note that 



where 



for the £ = 2 mode, there is an ill-defined term due to the 
factor {£ — 2) in the denominator. This is why we need a 
special treatment for the £ = 1 mode in the next section. 



4- For lower £ modes 

In the calculation of the second-order £ = 2 odd parity 
perturbation, we have the first-order £ = 1 mode con- 
tribution. In [38], this £ = 1 mode has been calculated 

under the Zerilli gauge, i.e., K\ m = h^ m = h[ e ^ m = 0. 



H$l(t,r) = — 1 ^ (r 3 ^f m (t) + M f m (t) 

oimV , J 3M(r-2M) 2 \ dt 2 ' ' v ' 1 w 



x6(r-R[t)), 
r d 
~(r- 2M) 2 dt 



f m (t)9(r-R(t)) 



H^(t,r) = {r _\ M)2 f m {t)e{r-R{t)), (33) 



f m (t)=SwnU(t) 



(R(t)~2M) 2 
R(t) 



Y^Qo, $(*)). (34) 



Here * denotes the complex conjugate. There is no con- 
tribution from the m — mode. 

Using the above first-order £ = 1 mode, we calculate 
the second-order source term from the coupling between 
this mode and the black hole's spin. Then the source 
term becomes finite at the horizon. In order to remove 
this finite term, we introduce a regularization function, 



(o,Z,2) 
2m 



(t,r)=* 



(o,Z,2),R 
2m 



SVl5(2-m) (2 + to) 
30 M r (r - 2 M) 



(t,r) 

f m {t)6{r-R{t)) , (35) 



and we solve the wave equation for the regularized func- 
tion ^E , 2m Z ' 2 ' ) ' R - Here, we note that the regularization 
function does not affect the waveform at infinity in our 
calculation. The regularized second-order source term is 
derived as 
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15 



(R(t)-2M) 2 U(t) ($(*))" 



R(t)M 



im(R(t)-2My$(t) R(t) (R(t)-2M) \ d 



(i) £7 (i) (i? (i) - 2 M) ($ (t)) 

R{t)M {R{t)fMU{t) ' (R(t)f MU{t) J dr 

™ (*) - 2 M ) ( 13 M - 3 i2 (t)) (* (t)) 3 U (t) 2 im (5 M — 2 i? (<)) (i? (t) - 2 M ) 2 $ (t) £7 (t) 



(R(t)) M 



(R (t)) M 



(i? (i) - 2 M) 2 i? (t) C7 (t) (- 12 M + 2 - m2i? (*) + 4 A (*)) (* (*)) R (*) (*) 



(R(t)) 5 M (R (t)) M 

im(llM - 2R(t))(R(t) -2M)6(i) (4 M - R (t)) R (t) 



(R (t)) MU (t) 



(R(t)) MU (t) 
I 



*(r - ii(t)) 



(36) 



We have only the local source contributions as the 
second-order source term from this mode. Using the 
following asymptotic behavior near the horizon, U(t) ~ 
(1 - 2M/R(t))- 1 , R(t) - (1 - 2M/R(t)), and $(t) - 
(1 — 2M / R(t)) , we find that the above source term van- 
ishes at the horizon in the integration of the wave equa- 
tion. 



(o) 



5. Symmetry in and ^ 



In this subsection, we use the notation 'pfcT 61 ^ = ' i 



and * 



(odd) 
im 



^imi which have the following relation in 



the first perturbative order: 



(even/odd) 



(-i) m Ut 



(cvcn/odd) 



(37) 



This is derived from a formula for the spherical harmon- 



Y e _ m (e,<p) = (-i) m Yi m (e,, 



(38) 



In the 0(a 1 ) calculation, we should have the same sym- 
metry because the metric perturbations become real. We 
can check this by using the explicit form of S^ en dd \ 



Gravitational waves 



In the above sections, we discussed the techniques to 



calculate the wave functions = 



Xl>(2) ^(o,l) 
im * £m> * im 



and The first-order wave functions and wave- 



,(o,Z,2) 

" Im 

forms at infinity arc simply related as 

y/(£ -!)£(£+!)(£ + 2) 



-ih x = 



2r 



where -QXtm denotes the spin-weighted spherical har- 
monics used in [42] . 

On the other hand, in order to discuss gravitational 
waveforms in the second perturbative order, we need to 
check the asymptotic behavior of the metric perturbation 
and the contributions from the first-order gauge trans- 
formation. First, we evaluate the asymptotic behavior 

of the tensor harmonics coefficients of G^J , because this 
information is used to construct the metric perturbation 
from the wave functions. For the odd parity-spin cou- 
pling part, we have the following behavior. 



x 0im 

4 2 2(o, S) 



B 



0(l/r 3 
0(l/r 3 
0(l/r 3 
0(l/r 4 
0(l/r 3 

and for the even parity-spin coupling part 



0(2) 
^■im 



(0,S) 



-0(l/r 3 ), 
^0(l/r 3 ), 
0(l/r 3 ), 
0(l/r 3 ), 
0(l/r 4 ) , (40) 





- 0(l/r 3 ) 




~ 0(l/r 3 ) , 




- 0(l/r 3 ) 




^0(l/r 2 ), 




- 0(l/r 2 ) 




0(l/r 4 ), 




- 0(l/r 3 ) 




^0(l/r 2 ), 






^0(l/r 3 ). (41) 



And the even parity-spin coupling part from the £ = 1 
even parity has a different behavior. 



Q®iE,S,[l = l]) ~ 0(l/r 3 ) 
Q<?l n (E ) S ) [e=l]) ~ 0(l/r 2 ) 



(42) 



(39) and V ( ^(E, S, [£ = 1]) = in the first-order Zerilli gauge. 
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From the above asymptotic behaviors, if we set the ob- 
server location to a large distance, we do not need to con- 
sider these tensor harmonics contributions because the 
contributions are at least 0(l/r) lower than the leading 
part. Note that the metric reconstruction in the second- 
order odd parity perturbation has been done from the 



Zerilli waveform ^ 



(o,Z,2) 



Next, we discuss the contributions from the first-order 
gauge transformation. Formally the following gauge 



transformation 
tion. 



is used in the second-order calcula- 



J RW 



V RW 



+ (x a ) + - \^ 2) ^ (x a ) + (x a ) 



r 



(43) 



where comma " ," in the index indicates the partial 
derivative with respect to the background Schwarzschild 
coordinates, and and £( 2 )^ are generators of the 

first and second-order gauge transformations, respec- 
tively. The subscripts RW and AF show the Regge- 



Wheeler gauge where we reconstruct the metric pertur- 
bation, and the asymptotic flat gauge where we obtain 
the gravitational waveforms, respectively. Then the met- 
ric perturbations change to 



L RW(iv 



-> h 



h 



(2) 

RWfiv 



(1) 

AFp,u 



r n 



,(2) 
l AF i iv 



l (2) 



r 



(44) 
(45) 



where C^i) denotes the Lie derivative. 

In this paper, second perturbative order means O(jia) 
where /i and a are small quantities. Since is 0(n), 
we ignore Cl^g^ and C^x)h^ Wllv with lvg Wllv ~ O(fi) 
On the other hand, there is a contribution 



in Eq. (45 I 



from £tm h 



(l,spin) 



The asymptotic behavior of this ten- 



sor harmonics coefficient becomes 



5H 0im ~ 0(l/r) , SH iem ~ 0(l/r) , 



(e) 

Olm 



0{r°) 



SH 2 em = , 5h, 

~ O(r ) , 8G lm ~ 0(l/r 2 ) , 
5K tm ~ 0(l/r 2 ) , Shotm ~ 0(r°) , 
Sh Um ~ O(r ) , Sh 2lm ~ 0(r°) . (46) 

For the I = 1 mode in the first perturbative order when 
we consider the gauge transformation to the center of 
mass coordinates, we have the same behaviors. These 
contributions to the second-order metric perturbation 
under the Regge- Wheeler gauge are also lower order by 
0(l/r) at least. 

Finally, to derive the waveforms in the SRWZ formal- 
ism, we may consider 



y/(e-me + w + 2) 



2r 



x $ 



lm 



(o) 



(47) 



where again $ tm = + *im- Note that for *im we 
have used a different wave functions for the first and sec- 
ond order odd parity calculations for the sake of simplic- 
ity of the final results. Using Eq. (31) and the above 

(2) . 

asymptotic behaviors of Q,\', we simply combine them 

as 



(o) 



= * 



(o,l) 
lm 



(o,Z,2) 
lm 



(48) 



7. Observer location effect 

In [IT], we saw that the observer location effect was 
not negligible on the waveforms. To compare the NR and 
perturbative waveforms, we directly use Eq. (47) because 



we can set the extraction radius of gravitational waves 
at a sufficiently distant location, for example, Rohs/M = 
1000. On the other hand, the NR waveforms are obtained 
from the NR -04 data 



tpi = h + — i h y 



(49) 



We should note that these are true only at i?obs — > oo. 

First, we discuss the asymptotic behavior of the (first- 
order) Regge- Wheeler-Zerilli functions. In general £ 
modes for both the even and odd parities, which we de- 
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note by and > are given by 



^vcn/odd) = HUt _ rl + ^±ll J dtHtm {t-r*) 



+0(r- 2 ). 



(50) 



We note that errors due to finite extraction radii, which 



arise from the integral term in Eq. (50), are larger for 



lower frequencies due to the 1 /u> factor obtained by inte- 
grating a function with frequency u). 

Next, we discuss the relation between Regge-Wheeler- 
Zcrilli functions and the mode function ipl m of the Weyl 
scalar. Here, although we can use the formula given in 
Eqs. (C.l) and (C.2) of [33], we use simpler formulae 
valid for the asymptotic behavior of the functions. If the 
NR Weyl scalar satisfies the Teukolsky equation in the 
Schwarzschild spacetime, we have 



r^ 4 m = H hn (t-r* 
+0(r- 2 ), 



-H trn {t - r ) 



2r 



(51) 



where the difference between Hg m and Hg m in Eq. ( 50 ) 
is only the numerical factor. 



Combining the above equation with Eq. ( 50 ) , we have 



ripj 



(even/odd) 



1 



(even/odd) 



-0(r" 



(52) 



which is independent of £ and parity modes. This equa- 
tion is consistent with the formula in [33]- Here, we have 
considered the correction for the RWZ functions. It is 
important, however, to calculate Hg m , the waveform at 
infinity, because the PN waveforms which are used to 
construct the hybrid waveform, do not have the finite 
observer location effects. 

Therefore, we consider the extrapolation of the NR ip^ 
from for example, Rq^ s /A1 — 100 to infinity by using 
Eq. Jsil: 



Him — 



em {£-!){£ + 2) 

-KObs %>A o 



0{R, 



Obs) ' 



(53) 



Again, the above formula is derived by assuming the 
Teukolsky equation in the Schwarzschild spacetime (a = 
0). Since we treat only the extrapolation from Robs/M = 
100 to infinity, we may use the wave (linear propagation) 
equation in the flat spacetime. Thus, the Teukolsky equa- 
tion with M —> is sufficient to discuss the extrapolation. 



This calculation gives the same result as Eq. ( 53 1 . Note 



that since the above formulation has been discussed by 
using the Weyl scalar in the Kinnersley tetrad, we need 
an extra factor as the explanation below Eq. ^ for that 
in another tetrad. 

Let us point out that full numerical methods using 
Cauchy-characteristic methods have been developed [45] . 
Also multipatch [35] and pseudospectral [TU] techniques 
allow extraction radii very far from the source. 



8. Numerical integration method 



Although we have used the combination of Eq. (27) 
for the even parity perturbation and integrate Eq. (28) 



in this paper, the basic equations are the four wave equa- 
tions, ([22]) and (p3l for the first perturbative order, and 



(25) and (|29fo for the second perturbative order. 



In order to integrate the resulting even and odd par- 
ity wave equations, we use the method described in [37]. 
This method is second-order accurate in the grid spacing 
(see [3T] for a fourth-order formalism) , but deals with the 
Dirac's delta source "exactly" or as accurately as needed. 

Even if we considered the metric ^ with first-order 
spin corrections to the Schwarzschild metric, the method 
of perturbations we used still propagates waves on the 
exact Schwarzschild background and lumps the spin cor- 
rections in a source term, as if they would be second- 
order perturbations. We hence apply the methods of 
[311 137] with an added smooth source to integrate the 
first-order in spin corrected RWZ wave equations. We 
proved second-order convergence of the extracted wave- 
forms and used spatial and time steps that produced er- 
rors well below those acceptable for full numerical evolu- 
tions. The runs typically take under a minute on a laptop 
and are very low in memory and resources requirements. 
We also note that these types of codes are amenable to 
implementation on accelerated hardware such as CPUs 
or Cell processors [4*5] , 



IV. ANALYSIS OF THE NUMERICAL VERSUS 
PERTURBATIVE RESULTS 

Here we directly compare the waveforms generated 
fully numerically with those computed by the perturba- 
tive (SRWZ) approach. Since our full numerical evolu- 
tions routinely extract the Weyl scalar ip4 at interme- 
diate radii, typically around R = 100M (a compromise 
between far enough of the sources and high enough lo- 
cal resolution), and the perturbative code evolves the 
Regge- Wheeler and Zerilli waveforms, we need to trans- 
late these different measurements of the waveform into a 
common radiation quantity. While analytic expressions 
already exists that relate them both [33J , such expressions 
involve second derivatives that lead to some numerical 
noise when building up tp^, for instance. The usual strain 
h also involves two integration constants that are hard to 
fix with accuracy [§1I32]- Hence, as a compromise, we use 
the news function, essentially dh/dt, which displays nicer 
smoothness properties for numerical comparisons. 

In Figs. [7j9]we superpose the waveforms obtained for 
the full numerical evolution of the q = 1/10 black-hole bi- 
nary case and the perturbative waveforms as computed 
by the integration of the wave equations ( p5| and (29) 
both, including the spin corrections (a/M = 0.26) or 
simply setting it to zero. We do these comparisons for 
the leading (£,m) — (2,2) mode and the next to leading 
(2,1) and (3,3) modes. Note that while (2,1) is an odd 
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FIG. 7: The real part of the (I = 2, m = 2) mode of dh/dt for 
the g = 1/10 case. The (black) solid, (red) dotted, and (blue) 
dashed curves show the NR, spin-off, and spin-on calculations, 
respectively. 




FIG. 8: The real part of the (£ = 2, m = 1) mode of dh/dt for 
the q — 1/10 case. The (black) solid, (red) dotted, and (blue) 
dashed curves show the NR, spin-off, and spin-on calculations, 
respectively. 
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parity mode (for a = 0) and comes from integration of the 
Regge- Wheeler equation (23), the other modes are even 



parity and hence obtained by integration of the Zerilli 
equation (22 1. In all cases we use the same "full numeri- 



cal" trajectory. When spin terms are switched on, there 
is a coupling of even and odd parity modes as shown in 
Eqs. p5| and p9l). 



We have computed the overlap functions, as defined 
in Ref. [9 , of these three sets of waveforms in order to 
quantify the phase agreement between them. This pro- 
vides some insight into the possibility of using these per- 
turbative waveforms to build up a bank of templates to 
support detection and analysis of gravitational wave ob- 
servatories such as LIGO and VIRGO. Table IIVI shows 
that the agreement between numerical and perturbative 
waveforms is very good in general for all three modes, and 
that including the spin dependence improves the match- 
ing to an excellent level. This improvement is based on 
the accurate description of the late time phase, as we will 
discuss next, and is independent of the particle's track. 
The orbital (inspiral) part of the waveforms are not so 
strongly dependent on the spin terms (for our simula- 
tions) and are correctly described by the nonspinning 
perturbations. It is interesting to note here that the ex- 
cellent phase agreement during the inspiral orbit might 
not be so surprising since the perturbative code uses the 
full numerical tracks (transformed into Schwarzschild co- 
ordinates); however, coordinates and gauges in full nu- 
merical evolutions are described in quite a different way 
than in (analytic) perturbative expressions and it is reas- 
suring to find such a good agreement in the final products 
of evolutions. 

In Figs. [10]|T2| we superpose the waveforms for the 
modes (2,2), (2,1), and (3,3) obtained from the full nu- 
merical evolution of the q = 1/15 case. We included full 
numerical, perturbative with spin (a/M = 0.189) and 



FIG. 9: The real part of the (£ = 3, m = 3) mode of dh/dt for 
the q — 1/10 case. The (black) solid, (red) dotted, and (blue) 
dashed curves show the NR, spin-off, and spin-on calculations, 
respectively. 




without spin corrections (a — 0) . We computed the over- 
lap functions, as defined in Ref. [3], for these three sets of 
waveforms and display the results in Table [Vj We observe 
again the generally very good agreement of the pertur- 
bative and full numerical waveforms. The agreement is 
still stronger when we include the spin dependence of the 
remnant black hole. 

In order to study in more detail the agreement of the 
numerical and perturbative waveforms we will proceed to 
decompose them into phase and amplitude (ip, A) with 
the usual formula 



ip = A exp(i(p) . 



(54) 



We display in Figs. pp5] the phases of the (2,2), 
(2,1) and (3,3) modes for the q = 1/10 case. Note the 
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TABLE IV: The overlap (matching) between the NR and perturbative dh/dt for the q = 1/10 case. The integration time is 
from t/M = 100 to 1220 and the definition of the matching is given in Eqs. (26) and (27) of [9]. 



Mode 


K(^ = 2,m = 


2) 


K(£ = 2,m = 


1) 


%t(£ = 3,m = 


3) 


Match (Spin OFF) 


0.980404 




0.968137 




0.927807 




Match (Spin ON) 


0.995055 




0.982173 




0.995347 




Mode 


Q(£ = 2,m = 


2) 


Q(£ = 2,m = 


1) 


Q(£ = 3,m = 


3) 


Match (Spin OFF) 


0.980379 




0.972727 




0.928151 




Match (Spin ON) 


0.995196 




0.982604 




0.995571 





TABLE V: The overlap (matching) between the NR and perturbative dh/dt for the q = 1/15 case, 
from t/M = 100 to 750, and the definition of the matching is given in Eqs. (26) and (27) of [9]- 



The integration time is 



Mode 


ft(£ = 2,m = 


2) 


"St(£ = 2,m = 


1) 


St(£ = 3,m = 


3) 


Spin OFF 


0.991297 




0.993986 




0.969254 




Spin ON 


0.996607 




0.997256 




0.995974 




Mode 


Q(£ = 2,m = 


2) 


Q(£ = 2,m = 


1) 


Q(£ = 3,m = 


3) 


Spin OFF 


0.991653 




0.996433 




0.968889 




Spin ON 


0.996780 




0.998178 




0.996218 





FIG. 10: The real part of the (I = 2, m = 2) mode of dh/dt for 
the q — 1/15 case. The (black) solid, (red) dotted, and (blue) 
dashed curves show the NR, spin-off, and spin-on calculations, 
respectively. 

0.03 ■ 



0.01 




FIG. 11: The real part of the (£ = 2, m = 1) mode of dh/dt for 
the q — 1/15 case. The (black) solid, (red) dotted, and (blue) 
dashed curves show the NR, spin-off, and spin-on calculations, 
respectively. 
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very good agreement between numerical and perturba- 
tive waveforms for the whole range of the simulation. All 
the agreements have been found with a single full numer- 
ical trajectory feeding the source terms of both the even 
and odd parity perturbative equations. The insets in the 
figures zoom in on the late time phases to display the ef- 
fect of the spin correction, which in all three modes shows 
improvements over the nonspinning background case. 

Figsures 16 18 show the phases of the (2,2), (2,1) and 
(3,3) modes for the q = 1/15 case. Again very good 
agreement is seen for the whole range of the full numerical 
simulation between perturbative and numerical results. 
The insets show that the spin correction, even if smaller 
than for the q = 1/10 case, still improves the late time 



phase, correctly capturing the quasinormal frequencies of 
the slowly rotating Kerr black hole (a/M = 0.189). 

We now turn to compare amplitudes of waveforms. Al- 
though for gravitational wave detection by the LIGO and 
VIRGO observatories the most important indicator is the 
phase, the amplitude agreement is particularly important 
in the modeling of the sources. Figure [19] directly com- 
pares the amplitudes of the q = 1/10 and q = 1/15 cases, 
shifted in time to agree at the peaks of their amplitudes. 
We then rescale the amplitudes of the q = 1/15 wave- 
form by the factor fi(q = 1/10)/ fi(q = 1/15) « 1.41 to 
verify a linear rescaling. We find that the rescaled am- 
plitude of the q — 1/15 wave is very close to the actual 
q = 1/10 amplitude showing that the systems are close 
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FIG. 12: The real part of the (£ = 3,m = 3) mode of dh/dt for 
the q — 1/15 case. The (black) solid, (red) dotted, and (blue) 
dashed curves show the NR, spin-off, and spin-on calculations, 
respectively. 
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FIG. 14: The phase evolution of the (£ = 2, m = 1) wave for 
the q — 1/10 case. The (black) solid, (red) dotted, and (blue) 
dashed curves show the NR, spin-off, and spin-on calculations, 
respectively. The inset shows the zoom-in for the quasinormal 
region. 




_i i i i i i i 

200 400 600 800 

t/M 



FIG. 13: The phase evolution of the (I = 2, m = 2) wave for 
the q — 1/10 case. The (black) solid, (red) dotted, and (blue) 
dashed curves show the NR, spin-off, and spin-on calculations, 
respectively. The inset shows the zoom-in for the quasinormal 
region. 
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FIG. 15: The phase evolution of the (i = 3, m = 3) wave for 
the q — 1/10 case. The (black) solid, (red) dotted, and (blue) 
dashed curves show the NR, spin-off, and spin-on calculations, 
respectively. The inset shows the zoom-in for the quasinormal 
region. 
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to behaving linearly at these mass ratios. 

In order to assess this last point in more detail, we 
compute the differences of the numerical and perturba- 
tive waveforms for each case, q = 1/10 and q = 1/15, and 
study how this "error" scales with q (or more precisely 
/j). We display the results of such computations in Figs. 
[20| and [21] for the cases of neglecting the spin of the fi- 
nal hole and that of taking it into account, respectively. 
The plots show that the inspiral phase scales like /i 2 as 
one would predict if the system would be completely lin- 
earized. While in the final merger region, near the peak 
of the amplitude, the rescaled differences display a de- 
pendence in fi between linear and quadratic, as if there 



are still nonlinearities present. One would expect this be- 
havior for values of q that are in the intermediate mass 
ratio regime, where the linear approximation is good but 
small nonlinear effects can still be observed. 



V. DISCUSSION 

In this paper we have described in detail the techniques 
used to compute gravitational waveforms with the per- 
turbative approach using full numerical trajectories in 
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FIG. 16: The phase evolution of the (I — 2, m = 2) wave for 
the q — 1/15 case. The (black) solid, (red) dotted, and (blue) 
dashed curves show the NR, spin-off, and spin-on calculations, 
respectively. The inset shows the zoom-in for the quasinormal 
region. 
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FIG. 18: The phase evolution of the (£ = 3, m — 3) wave for 
the q — 1/15 case. The (black) solid, (red) dotted, and (blue) 
dashed curves show the NR, spin-off, and spin-on calculations, 
respectively. The inset shows the zoom-in for the quasinormal 
region. 
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FIG. 17: The phase evolution of the (t = 2, m = 1) wave for 
the q — 1/15 case. The (black) solid, (red) dotted, and (blue) 
dashed curves show the NR, spin-off, and spin-on calculations, 
respectively. The inset shows the zoom-in for the quasinormal 
region. 
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FIG. 19: The amplitude of the (I = 2, m = 2) mode of the 
NR dh/dt for the q = 1/10 and 1/15 cases. The (black) 
thick solid and (red) solid curves show the q = 1/10 and 
1/15 amplitudes, respectively. The (red) dashed curve de- 
notes r](q = l/W)/r/(q = 1/15) ~ 1.41 times the q = 1/15 
amplitude. 
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the source terms of the perturbative wave equations. The 
program was successfully tested in the q = 1/10 case in 
Ref. [IT]. We have taken it further here studying larger 
initial separations for the full numerical evolutions of the 
q = 1/10 case, leading to simulations lasting for nearly 
eight orbits before the final plunge. We have also stud- 
ied the case q = 1/15, the smallest mass ratio so far 
in the literature, in order to assess quantitatively the q- 
dependence of the agreement of the full numerical and 
perturbative evolutions in the intermediate mass ratio 
regime. We have also included in our new computations 



the (linear dependence) spin of the final remnant in order 
to correctly reproduce the quasinormal ringing compo- 
nent of the full waveform at late times (after merger). 
The results are displayed in Tables |IV| and [V] and in 
Figs. |13||18| They show an apparent improvement in the 
matching (overlap) indices when the spin correction is 
taken into account compared to the vanishing spin case. 
In the Appendix we apply this linear-in-spin perturbation 
theory (SRWZ) to compute the corresponding quasinor- 
mal modes and compare the frequencies of these modes 
with those obtained for a Kerr black hole background 
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FIG. 20: The amplitude difference in the [i = 2, m = 2) mode 
between the NR and perturbative dh/dt for the spin-off cases. 
The (black) thick solid curve shows the q = 1/10 case. The 
(red) solid, dotted, and dashed curves show the amplitude 
differences for the q = 1/15 case rescaled by factors of 1, 1.41, 
and 1.41 2 , respectively. 
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FIG. 21: The amplitude difference in the (£ = 2, m = 2) mode 
between the NR and perturbative dh/dt for the spin-on cases. 
The (black) thick solid curve shows the q — 1/10 case. The 
(red) solid, dotted, and dashed curves show the amplitude 
differences for the q = 1/15 case rescaled by factors of 1, 1.41, 
and 1.41 2 , respectively. 




for all values of the spin parameter. We observe the re- 
sults in Figs. [29] and [3U] They show that SRWZ provide 
reliable predictions for a/M < 0.3, which justifies their 
use for the cases studied here where a/M — 0.26,0.19 
for q = 1/10,1/15 respectively. The generalization to 
arbitrary spins requires solving the Teukolsky equation 
instead of the RWZ ones [50]. Note that the relevant 
spin-effects on the waveform are due to the spin of the 
large black hole, while the effects of the spin of the small 
hole on radiation will tend to be negligible as q decreases. 



The use of numerical trajectories to describe the motion 
of the small hole in the field of the larger one already 
incorporates the spin dependence where the effects are 
stronger. 

After comparing the perturbative and full numerical 
waveforms and verifying the accuracy of the former, there 
remains the question of accurately modeling the trajec- 
tories for small q BHBs. We have stressed here an im- 
portant fact, that the trajectory dependence disappears 
from the perturbative formulation once the black holes 
merge, reducing the need of further full numerical sim- 
ulations with the resulting saving of computational re- 
sources. This savings is not negligible, because one not 
saves not only the (relatively short) time of evolution 
from merger to the end of the ringdown, but also the 
evolution time required to propagate the signal to an ob- 
server located far away from the sources. Typically, this 
should save over 500M of full numerical evolution. One 
can also predict the parameters of the final black hole by 
using formulae for the remnant parameters, as in |51U52) . 
found by empirical fitting. Still, the goal of our project is 
to be able to model, empirically, the BHBs inspiral tra- 
jectories as a function of q from a reasonably small num- 
ber of full numerical evolutions. In particular, numerical 
evolutions start from a finite, relatively close initial sep- 
aration of the holes. It is hence important to provide 
the large separation input from PN theory. While the 
full modeling of trajectories is beyond the scope of the 
current paper, here we discuss how this interface can be 
achieved for the current simulations of the q = 1/10 and 



q = 1/15 cases. The results are summarized in Figs. 22 
and [23] We have considered the full numerical and PN 
trajectories in the Schwarzschild coordinates, i.e., correct 
the full numerical tracks for the 1+log time slice and the 
PN ones for the quasi isotropic coordinates (ADM-TT 
gauge). In the q = 1/10 case, the full numerical evolu- 
tions essentially start from initial separations Ri « 9.5M 
in the Schwarzschild coordinates. We see a relatively 
smooth matching for the tracks and their first deriva- 



tive in (upper-left inset) Fig. 22 This would lead hence 
to smooth waveforms in the whole range of the evolu- 
tion, i.e., from as large initial (PN) separations as needed 
down to the ringdown. Note however, that in order to 
achieve this smooth matching of trajectories we had to 
make use of the resummed PN (RPN) evolutions (i.e. 
containing exactly the particle limit in the Schwarzschild 
background) . The RPN Hamiltonian used here is derived 
in the following. Based on the Hamiltonian formulation 
for the test particle given in [53] . the resummed part 
-ffsch is calculated by using the Schwarzschild metric in 
the isotropic coordinates. Then the RPN Hamiltonian is 
given by 



rrRPN _ tt 
n — -Hsch 



Hi 



PN 



H- 



2PN 



3PN 



(55) 



The finite mass effects iiipN, ^2PN and i?3PN in the 
above Hamiltonian arc introduced by the result of the 
standard 3PN Taylor Hamiltonian (TPN) and the 3.5PN 
radiation reaction effects on the equations of motion are 
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treated as in [S3]. In practice, i? n pN is obtained by the 
subtraction of the test particle limit from the Taylor PN 
Hamiltonian H n p^- The PN evolutions in the figures 
have been obtained from the quasicircular initial param- 
eter at R(t) ~ 50M. A good matching, at this initial 
separation, cannot be achieved with the TPN Hamilto- 
nian. Of course, at larger separations both PN expres- 
sions get closer to each other and a full numerical sim- 
ulation started at such large initial separations could be 
matched by Taylor PN expansion as well. 



FIG. 22: The radial trajectory R(t) obtained from the PN 
and NR evolutions for the q = 1/10 case in the Schwarzschild 
coordinates. The (black) solid, (red) dashed and (blue) dot- 
ted curves show the NR, resummed and PN Taylor ones, 
respectively. From the lower-right inset, we can choose the 
matching radius between the NR and resummed PN evolu- 
tions as R(t)/M = 9.35123. The upper- left inset is the zoom- 
in around the matching time. 
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This also suggest that, at even closer separations, as in 
the case of the numerical evolutions for q — 1/15 start- 
ing from Ri s» 8AM, not even the resummed PN leads 
to a very smooth matching of track. This is indeed the 
case displayed in Fig. [23] We may then conclude that, 
in order to simulate full inspirals of q ~ 1/10 matched 
to resummed PN, one needs to start the full numeri- 
cal simulations from initial separations Ri > 9M in the 
Schwarzschild coordinates, i.e., R^ 1 ^ > 8M in the quasi 
isotropic coordinates. Alternatively, one could seek to 
improve the resummed PN expansions with the effective- 
one-body (EOB) formalism [55] and its extension to in- 
corporate full numerical results (EOBNR) 56J . It is also 
relevant to cite here the works [5TH59"] that make pertur- 
bative evolutions of particle trajectories completely de- 
rived from PN expansions and used all the way down to 
merger without direct input from full numerical trajec- 
tories. 

If one indeed can extend those improved post- 
Newtonian treatments down to the ISCO in the particle 
limit, at R = 6M in the Schwarzschild coordinates, i.e. 
.R(iSCO) sa 4.95M in the isotropic coordinates, then one 



FIG. 23: The radial trajectory R(t) obtained from the PN 
and NR evolutions for the q = 1/15 case in the Schwarzschild 
coordinates. The (black) solid, (red) dashed and (blue) dot- 
ted curves show the NR, resummed and Taylor PN ones, 
respectively. From the lower-right inset, we can choose the 
matching radius between the NR and resummed PN evolu- 
tions as R(t)/M = 8.28796. The upper-left inset is the zoom- 
in around the matching time. 




can argue that the subsequent merger trajectory reaches 
a "universal" limit given by the geodesic motion of qua- 
sicircular orbits. In fact this seems to be the case for the 
tracks of the q — 1/10 and q — 1/15 simulations as dis- 
played in Fig. [3] One can argue that the very low level 
of radiation of those plunging orbits implies the univer- 
sal form of the track. This was also recently observed in 
[59] studying PN orbits. Notably, at the other extreme 
of the mass ratio range, i.e. for equal (and comparable) 
mass BHBs the strong gravitational emission taking place 
during the plunge erases any details of the preliminary 
evolution and one observes a universal waveform [60H63] 
To see the universal behavior of geodesies inside the 
ISCO for quasicircular inspirals, we use the orbits with 
imaginary eccentricities for timelike geodesies in the 
Schwarzschild spacetime as given on page 111 of |64j . 
The initial part of these orbits can be considered the con- 
tinuation of the inspiral trajectories through the ISCO. 
These geodesies have the following form near the horizon: 



*(JZ) 



3 fi-AY 

6M J 




(56) 



where the imaginary eccentricity (ie) is a small quantity, 
and R <6M. 
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The initial velocity at R(t) 
liven by 



Rq is approximately 



dR(t) 
dt 

d$(t) 
dt 



V6 



= — 7t- e 



V6 



6 

V6 

36M + 24M 

V6 



6M J ' 

1-*) 
6M J 



+ 



288M 



15 1 



-RoY 
6M I 



(57) 



which allows us to match to full numerical trajectories 
and then use the geodesic expressions to smoothly sup- 
press the local source terms when the particle approaches 
the Schwarzschild horizon (see Sec. 



Ill A 4 ) 



In Fig. [24} we plot the phase evolution in terms of the 
orbital radius. As a fiducial starting point, just inside 
the Schwarzschild ISCO, we take the self-force corrected 
ISCO radius 



R = 6 M- 3.269 fi, 



(58) 



as discussed in [65j . Although we see some differences 
in the initial part of the orbits, the trajectories reach a 
universal limit approaching the horizon. 



FIG. 24: The orbit with imaginary eccentricities discussed 
in 64 . The thick and thin curves show the q = 1/10 and 
q = 1/15 cases, respectively, 
various eccentricities. 

1 r 



Here we show the orbits with 
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Appendix A: Analysis of the wave equations 

The following is useful for the analytic discussion, espe- 
cially the behavior of the source term around the horizon. 
And ths also gives a stable evolution in the numerical cal- 
culation, because nonvanishing contributions at the hori- 
zon in the source terms are canceled out analytically. 

Here, we discuss the wave functions as 

*A»(*,r) - ^(t,r)6(R(t)-r) 



(out) 



(step) 



(t,r) 



(out; 
em 



(t,r)9(r-R(t)). 



(Al) 



where ^>t m denotes the even or odd parity wave function. 



The functions *ftf 5 



and *^ ep ^ are the homoge- 



neous solutions of the Regge-Wheclcr-Zerilli equations. 
From these definition, we have 



-*f m ep) (t,r)6(r-R(t)). (A2) 



Therefore, for example, the time derivative of the above 
wave function is written as 



(step) 



(t,r) 0(r 



(rf»p) (ijr) ^RW 5(r 



'd t *t ep \t,r))8(r-R(t)) 



dt 



R(t)) 
-R(t)) 



r - R{t)) . (A3) 



To find the quantities of the waveforms at the particle 
location, i.e., ^^ cp ^(t, R(t)), we use 



8 vI/W 



(t,r) 



(z,i) 

till 



(t,r) 



16 V^nir 2 (r - 2M) (i) 
'£(£+1) (r£ 2 + rl - 2 r + 6 M) Um 



16v / 27rzr(r-2M) {1) 



(i -!)(£ + 2) ^W+T) 



e&(*»r), 



(A4) 



where each wave function in the left-hand and right- 
hand side of the above equations behaves as a step func- 
tion at the particle's location because of the first-order 
Regge-Whcclcr-Zcrilli waveforms. Therefore, substitut- 



ing Eq. (IA3J into and d t ^ lni 



(o.l) 



we obtain the an- 
alytic expression of ^f m ep \t, R(t)) from the coefficients 
of the Dirac's delta function. 
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1. Analysis of the even parity wave equation 

We have introduced a new function for the even parity 
calculation to the SRWZ formalism, 

*t m (t,r) = *£(t,r) + *%l(t 1 r) . (A5) 

The gravitational waveform with the spin effect is ob- 
tained directly from ff m . Therefore, we discuss the 



wave equation for $| m in the following. Here, we cre- 
ate our numerical code for the perturbative calculation 
based on (3Tj . It is important to distinguish the cell that 
the particle does cross from the other cells. 

For the cell that the particle does not cross, we use 
the following homogeneous equation, i.e., can read the 
following equation from the step function part, which 
does not include the local source term: 



d 2 T , x (r-2Af) 2 <9 2 , . „(r-2M)M9 T . . 
- W ^lm {t, r) + r2 ^ tm (t, r) + 2 ^ -3-^— —9 im (t, r) 

- (r - 2 M) (4 r 3 £ - £ 4 r 3 + 3 l b r 3 - 7 £ 3 r 3 + £ 6 r 3 + 12 £ 3 r 2 M - 24 r 2 Ml - 18 r 2 M£ 2 + 24 r 2 M 

+6£ 4 r 2 M - 72 rM 2 + 36£ 2 rM 2 + 36£rM 2 + 72 A/ 3 ) [t, r) / Ur£ 2 + £r ~ 2 r + 6 M) 2 r 4 

-4 iS m (4 r 3 £ 7 + 144 A/ 3 ^ 2 + 16 r 3 £ - 24 r 3 + 18 A«V + 144 Af 3 £ + r 3 £ 8 - 216 rM 2 - 66 fV 2 Af 
-48r 2 A« + 144 M 3 + 36^ 2 rAf 2 + 22 r 3 £ 2 + 120 r 2 M + 6fr 3 - 11 £ 4 r 3 + 54Mr 2 £ 5 
+72 M 2 £ 4 r + 12 £ 4 r 2 M + lUM 2 r£ 3 - 90 r 2 M£ 2 - 36£rM 2 - U£ 5 r 3 ) J^m (*, r) 



I [r 3 {£+!)£ (r£ 2 +£r-2r + 6M) : 



24iS (I + 2) (I- l)m(r-2A/) 2 d 2 ^ 
r 2 (r£ 2 + £r - 2 r + 6 M f £ {£ + 1) 9tcV 



12rA/£ 2 



12 5 w if„ m )M] m \ (r - 2 M ) {£ - 2) (£ 5 r 2 + 2 r 2 £ 4 + 4 r£ 3 M - 4 r 2 £ 2 + 
Y (2^—1) (2^ + 1) v ' v ; v 

-5r 2 £+ 12rM£ + 6r 2 - 28 rM + 36M 2 ) *^ lm (f,r) / [r 5 ^ (r£ 2 + £r - 2 r + 6 M)' 

12 5 ^ ^^mXi^+S)^ (r ~ 2 M) {i + 3) ^ V + 3 + 4 ^ M + 2 + 2 
+32 rM - 8 r 2 - 36 A/ 2 ) 4^ m (t, r) / (£ + 1) r 5 (rf 2 + £r - 2 r + 6 M)' 



(A6) 



And then, we need the following local source term 
which is added to the right hand side of the above equa- 
tion, for the cell that the particle does cross: 



where the first-order source term S, 



,i,L) 



is the same 



as S t 
f3U as 



(even,l) 
Im 



in Section 



IIIB1 



and given in Eq. (A. 5) of 



£,(evcn.L) ^,(cvcn,l,L) 



-,(cvcn,2,L) 



(A7) 
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32tt M (R(t)-2M) (2 M — R(t) — R (t) R (tj) (2 M - R(t) + R(t) R(t)) U (t) j 

V L\ / —5(r - R(t)) 

£(£+1) (R(t)) 2 (R(t)£ 2 + R(t)£-2R(t) + 6M) dr 

^ / 32m 2 (fl(t)-2M)[/(t)(4>(t)) 2 64imR(t)(R(t)-2M)U(t)$(t) 
+ £(£+l)\ (£-l)(£ + 2) R(t)£ 2 + R(t)£~2R(t) + 6M 



16 (R (t) — 2 M) U (t) 

' 8M + 10M£ 2 + 10 M£- 3 R(t)£ 2 



(R (t) f 2 + R (t) £ - 2 R (t) + 6 M) (£ - 1) (t + 2) 

i6U(t) (ii(t)) 2 

+2 R(t)£ 3 +AR (t) +R(t)£ i -4R (t) £) + ^ '- ? 

; R{t){R{t)£ 2 +R{t)£-2R{t) + QM) 2 

x (-2 (R (t)f £ - (R (t)) 2 £ 2 + 2 £ 3 (R (t)) 2 + £ 4 (R [t)f + 12R (t) £ 2 M +12R (t) £M + 12 M 2 ) 
16 (R (t) - 2 M) U (t) (Biy m2 + UR ^ fM _ 2AR(t) M + 12R (t) £M 



(R (t)) 6 (R (t) £ 2 +R(t)£-2R(t) + 6 M) 
-2 (R (t)f £ - (R (t)) 2 £ 2 + 2£ 3 (R (t)) 2 + £ 4 (R (i)) 2 ) \s(r-R (tj) 



F/ m (6o,$(i)) . (A8) 
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The second-order local source term <S^ en ' 2 ' L ^ has the following expression: 



-,(even,2,L) 



192imS*7r M [/(t) (R (t) - 2 M) R (t) (£ 2 + £ - 2 to 2 ) ($(t)Y A 

, ^ '-Y* (Go, $(*)) f-<f(r - i2(t)) 

(£ + 2)(£-l)(£ + l) 2 £ 2 (i?(i)£ 2 +i?(i)£-2i?(i)+6M) 



+ 



-24irrag (£ + 2) (£ - 1) (R (t) - 2 M) 2 i^(ste P ) 
(i?(t)) 2 (fl(i)£ 2 + fl(i)^-2i?(t)+6M) 2 ^ + l) & ^ m 



(t,r) 



r=R(t) 



+imSTTnY; m (e ,$(i)) (-384 [l 2 +£-2m 2 ) (30 M 2 + 6 i? (i) £ 2 M + 6 R (t) IM 



-21 R{t)M-2l{R{t)) 2 -2£ 2 (R(t)) 2 + 4(R(t)) 2 )U(t) ($ (t)J £(t) 

/ [i? (t) (£ + 2) (£ - 1) (£ + l) 2 £ 2 (i? (t) ^ 2 + R (t) I - 2 R (t) + 6 M) 2 ] 

-128 (R (t) — 2 M) U (t) R (t) ((R (t)) 2 £ i + 2{R (t)) 2 £ 3 - 6£ 2 (R {t)f + 18R (t) £ 2 M 

-7£{R {t)f + 18 R (t) £M + 10 (R (t)) 2 - 36 R (t) M + 36 M 2 ) 

I [£ (R (t)) 3 (R (t) £ 2 + R (t) £ - 2 R (t) + 6 M) 3 (£ + 1)] 



+ - 



192 i m U (t) (£ 2 +£-2 to 2 ) (R (t) — 2M) \ § (t) 

(£ + 2) (I- 1) (t + l) 2 £ 2 (R (t) £ 2 + R (t) £ - 2 R (t) + 6 M) 
+128 imU(t) (R (t) - 2 M) 2 (72 M 2 + 30 R (t) £ 2 M + 30R (t) £M - 60 R (t) M 
-10£(R (t)) 2 + (R (t)) 2 £ 4 + 2(R (t)) 2 e + 16 (R (t)) 2 -9£ 2 (R (t)) 2 )<S>(t) 



I [{£ + l) 2 £ 2 (R (t) £ 2 + R(t)£-2R(t) + 6 My (R (t)Y 



5(r-R(t)) 



256 7r n S , 



{£ - to) {£ + to) U (t) $ (t) (R (t) — 2 M) 2 (—2 R(t) + R(t)£ + R (t) £ 2 + 3 M ) 



(2^-1) (2^+1) 



^ 2 + 1) (i? (i? (t) £ 2 + i? (i) £ - 2 i? (t) + 6 M) 



-256 7T /i 5. 



' + m + 1) {£ - to + 1) U (t) $ (t) (R (t) - 2 M) (-2 R(t) + R(t)£ + R (t) £ 2 + 3M) 



(2£+l)(2£ + 'S) 
xdgY; +lm (Q ,m)))s(r-R(t)) . 



£(£ + l) 2 (R (t)f (R (t) £ 2 + R (t) £ - 2 R (t) + 6 M) 2 



(A9) 



Here we have used the analytic expression of the wave 
function ^^ p \t, R(t)) at the particle's location, and 
the instantaneous geodesic approximation for the second- 
order source term. 

It is noted that there are remaining terms at the hori- 
zon in the integration of the wave equation. These arise 
from the transformation of the original wave equation to 
the above equation. Since the source term for the origi- 
nal wave equation docs not include any remaining term 
at the horizon, these remaining terms cancel out with the 
derivatives of the wave function, i.e., d^ v) /dt. 



2. Analysis of the odd parity wave equation 



In the calculation for the odd parity perturbation of the 
SRWZ formalism, we have treated the first- and second- 
order perturbations separately. The first-order (local) 
' v " L c ' 1 which has been sim- 



source term, 5. 



fm 



plificd with the geodesic equation, is given as 
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-,(odd,l,L) 



32 7T n 



£(£+!)(£ -l)(£ + 2) 



U(t) (R(t)-2M) {R{t)f 



(R (t) -2M)$(t) 



+ ((2R(t) — 7 M) R (t) U (t) ($(*)) -imR(t)R(t)U(t) 



[/ (t) 7 dr 

(i? (t) - 5 M) * (t) 
R(t)U(t) 



5(r-R(t)) 



2 (i?ft)-2M) 2 <7ft) 
(i?(t)) 2 



<f(r-fl(t)) 



a e r; m (e ,$W) . 



(A10) 



Next, we focus on the second-order wave equation. 
^±i TO an d have already been derived in the first- 

order calculation. For the cell that the particle does not 



cross, we may consider only the homogeneous part of the 
wave equation, 



~A°^ 2) M + ^<f' 2) (*>') - ^° dd) (0*t Z ' 2) (*,r) 



/ (2r£ 3 -5r£ 2 + 18r -6r£ + r£ 4 + 9M£ 2 -42M + 9£M) (r - 2 M) {ol) 



I / (l-m)(l + m) / (r-2Mr^-l) g 2 (1) 



3(r -2M)(£-l) (tr 2 - 2r 2 £ 3 ~r 2 £ 2 + 2r 2 £ + 6r£ 2 M -6r£M -l2rM + 24M 2 ) d {1) 
2 r 6 ( r £ 2 -r£-2r + 6M) di^ l - lm 



±1 l (£ + m, + l)(£-m+l) f (r-2M) 2 (£ + 2) d 2 (1) 



(£+l)(£ + 2)\/ (2£ + l)(2^ + 3) V r5 

2 f 3 + llr¥ 2 + 6r 2 £ + 6rf 2 M + 18.. , . , ,, 

r6 (r£ 2 + 3rf + 6M) ^+im(*> r )J- ( AU ) 



3 (r - 2 M) + 2) (£ A r 2 + 6 r 2 ^ 3 + f 1 r 2 £ 2 + 6r 2 £ + 6 r£ 2 M + 18 r£M + 24 M 2 ) d >T ( i) 



The second-order local source terms, S\ m ' ' ' which as 
we need for the cell that the particle does cross is written 



2G 



S, 



(odd,Z,2,L) 
Im 



V (R(t)f dr 

32vr/i (1 + 3) (£-2)U (t) $ (i) (R(t)-2M) 
(i? (£)) 5 (£+l) 2 ^ 2 + 2) 



r=fi(t) 
2 



-9 e Y/ ro (e Q ,<I>(t))U(r 



45 



- m) (£ + m) 
~^-l)V {2l-l){2l+\) 



I2nfi (£ + 2) (£ 2 -£-2m 2 ) ($ (t)) 2 (-i? (i) + 2 M ) f7 (t) i?(t) 



(fi(t)r(^-2)i 



r=R(t) 



tt M (£ + 2) / 12* m (i2 (t) - 2 M) C/ (f) (£ 2 - * - 2 m 2 ) ($ (i)) ' 



12 (5 R (i) - 14 M) R (t) U (t) (£ 2 -1-2 m 2 ) ($ (t)) ' 

(R(t)f(e-z) 

48 (R (t) - 2 Mf (£+l)U (t) R (t) {1-1)1 
(R (t)) 5 (R (t) £ 2 - R (t) £ - 2 R (t) + 6 M) 

o -£- 1 



96 i m (R (t) - 2 Mf_ (£+l)U (t) $ (t) 
(R (t)f (R (t) i 2 - R (i) I - 2 R (t) + 6 M) 



Y;_ lm (Q ,<S>(t)))5(r-R(t)) 



(A12) 



r 



where [£ -f-> — £ — 1] refers to an additional te rm ob tained 
by replacing £ with — £— 1 in all terms in Eq. ( A12 ) start- 
ing from jjj^p 



and subsequently replacing 'S^J^om 



with #^f„ p) and F^_ 2m with Y e * +lm . The above source 
term is added to the right hand side of the homogeneous 
part of the wave equation. It is found that the right hand 
side of the equation vanishes at the horizon. Here, the 
instantaneous geodesic approximation has also been used 
in the above equation. 



dominant and couplings with the other modes can be ig- 
nored. Also, the perturbed Regge-Wheeler-Zerilli equa- 
tions with the spin effect do not have the local source 
terms, i.e., we consider the homogeneous equation. 



3. Analysis of quasinormal modes 



In the SRWZ formalism, we discuss a special case 
where the £ — 2, m = ±2 even or odd parity mode is 



For the even parity part, we use the same equation as 
Eq. ( A6 1 without the other mode coupling, 



dt 2 



dr 



6 



(r - 2 M) (4 r 3 + 4 r 2 M + 6 rM 2 + 3 M 3 ) 



r i (2r + 3M) 2 



*2±2(t,r) 



, 8iS [r-2Mf d 2 T . , 8iS (6i 
^ 2 (2, + 3M)^ W ' r)T — 



46 r 2 M + 45 rM 2 + 21 M 3 ) d 
r 3 (2r + 3M) 3 dt 



2±2 



(t,r)=0. (A13) 



For the odd parity part, we use a different equation a simple equation by using only the Cunningham et al. 



from Eq. (|All|). This is because if we ignore the other waveform ^ { 2 % (or the Zerilli waveform * 2 ±1 } ). The 



mode coupling and the local source term, we can derive £ = 2, m = ±2 odd parity wave equation with the spin 
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effect becomes 



d 2 

2iS d 2 



6 



(r - M) (r - 2 M) 



± 



r 2 drdi 2±2 



2iS (7r 2 - 17rM + 8M 2 ) a 



(r - 2 M) r 4 
I 



0t 2±2 



(A14) 



where we have introduced \1/ 

2±2 (*) r )- 



(o,2), 



(o) 
2±2 



(o,l) 
2±2 



(*.r) 



We treat the above equations in the frequency domain, 



d 2 (r - 2 M)(4r 3 + 4r 2 M + 6rM 2 + 3 M 3 ) 



±- 



dr* 

iSuj(r-2M) 2 d 
r 2 (2r + 3M) 2 dr 



r 4 (2r + 3M) 5 



*2± 2 (^;t') 



*2±2 (w;r) =F 



; S u> (6 r 3 + 46 r 2 M + 45 r M 2 + 21 M 3 ) 



r 3 (2 r + 3 M ) 



dr* 



6 



(r - M) (r — 2 M) 



± 



25w d (o ) 
r 2 dr 2±2 



( W ;r)T 



25w (7r 2 - 17rA/ + 8Af 2 ) (o 



(r - 2 M) r 4 



*w a (w;r) = 0. 



#2± 2 (a;;r) =0, (A15) 



(A16) 



For the nonspinning (5 = 0) case of the above equa- 
tions, we have already known the transformation between 
the Regge- Wheeler and Zerilli function. This is known 
as the Chandrasekhar transformation 64], given by 

^ (tr) f 6 + 9 M2 (^~2M) \ (1) 

+ 3 m(i-2^) J^&ftr). (A17) 

Using these transformation, for example, we may solve 
only the Regge- Wheeler equation to obtain the quasinor- 
mal frequency. 

In order to discuss a similar treatment up to 0(a 1 ) 
(a = S/M), first we consider the following transforma- 



tion: 

*2±2(^;r) = exp fi^T^f) *2±2 (w; r) , 
tt<°2 2 ( W ; r) = exp U-^m) ("5 r ) > ( A18 ) 

where these transformations are consistent in the 0(a 1 ). 
Since we treat the wave functions only up to 0(a 1 ), we 
may choose another transformation here. From the above 
transformations, we have the simple differential equa- 
tions which are similar to the Regge- Wheeler and Zer- 
illi equations. The difference from the original Regge- 
Wheeler and Zerilli equations arises in the potential 
terms. 



d 2 (r - 2 M)(4 r 3 + 4 r 2 M + 6 rM 2 + 3 M 3 ) 
+ dr^~ ~ r 4 (2r + 3Af) 2 

8 S u) (4 r 3 + 56 r 2 M + 36 rM 2 + 15 M 3 ) - 



*2±2(w; 



r 3 (2r + 3M) d 

d 2 (r - M) (r - 2 M) 
— 6 



dr* 



* 2± 2 (w;r)=0, 

^;^ 45 - (3 r 2M) ^2(^;^o. 



(A19) 
(A20) 
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From these equations, we find the " Chandrasekhar" transformation as 



M 2 {r-2M) SMuj (45 M 2 - 48r 
r 2 (2r + 3M) T r 2 {2r + 3M) 



2 \ 



*2 2 (o.;r) = 6 + 9 " , - ' ~> T - — " " , ' 2 : ' | *2±2 for) 



The differential equations for the even and odd parity 
perturbation become the same form by using the above 
transformation. 

Next, we consider quasinormal modes derived from 



Eq. ( A20 ) . A recent review for quasinormal modes is 
given in |66) . Here, we should note that if we use 



Eq. (A18l to obtain the simple equation in Eq. (A20l, 



these change the boundary behaviors near the horizon 
and at infinity. Therefore, although the expression is 
same in the 0(a 1 ) expansion, we should consider to do 
another transformation: 



(o) 
2±2 



(w;r) 



1 



2M 



In 1 ± 



Slot 



{r-2Mf 



x*<°) 2 ( W ;r) . 



(A22) 



This does not change the boundary behaviors. 

In order to calculate the quasinormal frequencies, we 
use the Leaver's method [ST]. As boundary conditions, 
the wave function ^±2 nas the following behaviors: 



(o) 
2±2 

(o) 
2±2 



(p;r) — > r p e pr for r — » oo , 
(p;r) -> (r- l) p+lx for r ->• 1 



(A23) 



where we have considered 2M — I and p = —iu which 
are the same notation as [67]. Here, \ is defined by the 
nondimensional spin parameter \ — S/M 2 . Then a solu- 
tion of Eq. ( A20 ) can be written in the form of 



(o) 
2±2 



r -P e -p{r-l) ( r _l r - r 



x 22 a „ 



r - I 
r 



P+i \ r -{P+iX) 

(A24) 



We obtain the recurrence relation for a n in the above 
equation, 



a a\ + /3 a = , 

and for n > 1, 

a n a n +i + P n a n + 7„ a n _i 



(A25) 
(A26) 



where 

Pn 



(2 + 2 n) p + 2 i (n + 1 ) X + (n + 1) 



; p 2 + (-4 - 8 n - 7 ix) p - 2 i (2 n + 1) x 



-3- 2n 2 - 2n, 

4 p 2 + (4 n + 5 «x) /° + 2 i x n 

+ (n - 2) (n + 2) . 



(A27) 



When we set x = 0, the above equations reduce to Eq. (8) 
in [67]. 

In Fig. [25] we show the result for the quasinormal fre- 
quencies, u> around % = 0. As a reference, we also plot 
the values given in Table II of [55]. Figures 26 and 27 



show the real and imaginary parts of the quasinormal fre- 
quencies around x = 0, respectively. The Figures 28 29 



and [30] show the result for -0.9 <x< 0.9. In Table |VT 
we show the numerical values and the relative errors for 
the real and imaginary parts of p defined by 



ErrsR 



9E(p„) - Sft(p) 
R(p) 



Errc, = 



9f(p) 



(A28) 



where p a and p represent our result and that of [68j . 
respectively. We plot the above errors in Figs. [3X1 and 
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[32] and we zoom in the region —0.5 < x ^ 0-5 in Fig. 
which shows the absolute values of the relative error. 



FIG. 25: The quasinormal frequencies, oj around x = 0. We 
have used the same expression as |67| . The (red) circles show 
our result, and the + marks denote the values given in Table 

ii of m. 
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TABLE VI: The quasinormal frequencies in terms of p = 



The m = — 2 mode can be considered as the m = 2 mode with 



the inverse spin signature, 
calculation (see \ — 0-0). 



Here we set avr = in the recurrence relation of Eq. ( A26 1 . This creates the numerical error in our 
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m = 2 (Th 
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paper) 


m = 2 ([68 ) 


Errs 


Err a 


-0.9 


-0 


173072 - 





581783i 
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176562 - 
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-0.019766 
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-0 


177024 - 
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014255 
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011023 
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643379 i 
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648614* 


-0.006947 
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008071 
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680440 * 
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682666* 


-0.002455 


o 


003260 


-0.2 
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178181 - 





701019* 
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FIG. 26: The real part of the quasinormal frequencies, uj 
around \ — 0. The horizontal axis denotes the nondimen- 
sional spin parameter, \ — S/M 2 . The (red) circles show our 
result and the + marks denote the values given in Table II of 
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FIG. 27: The (minus) imaginary part of the quasinormal fre- 
quencies, uj around \ — 0. The horizontal axis denotes the 
nondimensional spin parameter, x = S/M 2 . The (red) circles 
show our result and the + marks denote the values given in 
Table II of 1551. 



0.1782 Q 



O 



O 



O 



O 



o 



o 



O 



O 



o 



o 



o 



o 



o 



o 



o 



o 



o 



6- 



0.05 0.1 

x 



FIG. 29: The real part of the quasinormal frequencies, uj for 
—0.9 < \ < 0.9. The horizontal axis denotes the nondimen- 
sional spin parameter, x = S/M 2 . The (red) circles show our 
result and the + marks denote the values given in Table II of 
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FIG. 28: The quasinormal frequencies, u> for —0.9 < x < 0.9. 
We have used the same expression of [67] • The (red) circles 
show our result and the + marks denote the values given in 
Table II of [68], 
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FIG. 30: The (minus) imaginary part of the quasinormal fre- 
quencies, u) for —0.9 < x < 0.9. The horizontal axis denotes 
the nondimensional spin parameter, x = S/M 2 . The (red) 
circles show our result and the + marks denote the values 
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FIG. 33: The absolute value of the relative errors in quasi- 
normal frequencies in the region —0.5 < x < 0.5. The (red) 
circles and (blue) boxes show those of the real and imaginary 
parts of the frequency, i.e., |Errg| and |Errsu|, respectively. 
The horizontal axis denotes the nondimensional spin param- 
eter, x = S/M 2 . 
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FIG. 32: The error in the (minus) imaginary part of the quasi- 
normal frequencies, i.e., Errjf . The horizontal axis denotes the 
nondimensional spin parameter, x = S/M 2 . 
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